#include "toroidal.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

void field::load_fromfile(const char* str)
{
    int i = 0;
    int j = 0;
    FILE *fp = NULL;
    char *line,*record;
    char buffer[20450];
    if((fp = fopen(str, "r")) != NULL)
    {
        //printf("file open success:%s\n",str);
        while ((line = fgets(buffer, sizeof(buffer), fp))!=NULL)
        {
            record = strtok(line, ",");
            while (record != NULL)  
            {
                this->data[i*Ny+j++] = atof(record);
                record = strtok(NULL, ",");
            }
            i++; j = 0;
        }
        fclose(fp);
        fp = NULL;
        //printf("file close success:%s\n",str);
    }
}

void field::output(const char* str)
{
    double* pointer = this->data;
    FILE* fp = fopen(str,"w");
    fprintf(fp,"R,Z,val\n");
    for(int i = 0; i < Nx*Ny; i++)
    {
        double R = i/Ny * dr + Rmin;
        double Z = i%Ny *dz + Zmin;
        fprintf(fp,"%f,%f,%f\n",R,Z,data[i]);
    }
    fclose(fp);
}

void field::vtkoutput(const char* str)
{
    double* pointer = this->data;
    FILE* fp = fopen(str,"w");
    fprintf(fp,"# vtk DataFile Version 2.0\nVolume example\nASCII\nDATASET STRUCTURED_POINTS");
    fprintf(fp,"\nDIMENSIONS %d %d 1\nASPECT_RATIO %f %f 1\nORIGIN %.2f %.2f 0\n",
    Nx, Ny, dr, dz, Rmin, Zmin);
    fprintf(fp,"POINT_DATA %d\nSCALARS volume_scalars float 1\nLOOKUP_TABLE default\n",
     Nx*Ny);

    for(int i = 0; i < Nx*Ny; i++)
    {
        int iy = i%Nx;
        int ix = i/Nx;

        fprintf(fp," %1.3f",data[iy*Ny+ix]);
    }
    fclose(fp);
}

void vfield::vtkoutput(const char* str)
{
    double* ptr = this->fr;
    double* ptz = this->fz;
    FILE* fp = fopen(str,"w");
    fprintf(fp,"# vtk DataFile Version 2.0\nVolume example\nASCII\nDATASET STRUCTURED_POINTS");
    fprintf(fp,"\nDIMENSIONS %d %d 1\nASPECT_RATIO %f %f 1\nORIGIN %.2f %.2f 0\n",
    Nx, Ny, dr, dz, Rmin, Zmin);
    fprintf(fp,"POINT_DATA %d\nSCALARS fieldr float 1\nLOOKUP_TABLE default\n",
     Nx*Ny);

    for(int i = 0; i < Nx*Ny; i++)
    {
        int iy = i%Nx;
        int ix = i/Nx;

        fprintf(fp," %1.3f",ptr[iy*Ny+ix]);
    }
    fprintf(fp,"\n\nSCALARS fieldz float 1\nLOOKUP_TABLE default\n",
     Nx*Ny);

    for(int i = 0; i < Nx*Ny; i++)
    {
        int iy = i%Nx;
        int ix = i/Nx;

        fprintf(fp," %1.3f",ptz[iy*Ny+ix]);
    }
    fclose(fp);
}


void fiq::output(const char* str)
{
    double* poi = this->f_inphase;
    double* poq = this->f_quadrature;
    FILE* fp = fopen(str,"w");
    fprintf(fp,"X,Y,Z,val\n");
    for(int i = 0; i < Nx*Ny; i++)
    for (double k = 0; k < 2*pi; k += 2*pi/32.0)
    {
        double R = i/Ny * dr + Rmin;
        double Z = i%Ny *dz + Zmin;
        if(checkboundary({R,0.0,Z}))
        fprintf(fp,"%f,%f,%f,%f\n",R*cos(k),R*sin(k),Z,
            poi[i]*cos(n*k)+poq[i]*sin(n*k));
    }
    fclose(fp);
}
void trace_output(const char* str, cylv* pts, int n)
{   
    FILE* fp = fopen(str,"w");
    fprintf(fp,"R,Z\n");
    for(int i = 0; i < n; i++)
    {
        vec a = cyl2xyz(pts[i]);
        fprintf(fp,"%f, %f\n", sqrt(a.x*a.x+a.y*a.y),a.z);
    }
    fclose(fp);
}

void trace_vtkoutput(const char* str, cylv* pts, int n)
{   
    FILE* fp = fopen(str,"w");
    fprintf(fp,"# vtk DataFile Version 2.0\n");
    fprintf(fp,"Cube example\n");
    fprintf(fp,"ASCII\n");
    fprintf(fp,"DATASET UNSTRUCTURED_GRID\n");
    fprintf(fp,"POINTS %d float\n",n);
    for(int i = 0; i < n; i++)
    {
        vec a = cyl2xyz(pts[i]);
        fprintf(fp,"%1.3f %1.3f %1.3f ", a.x, a.y,a.z);
        if(i % 13 == 12)
            fprintf(fp,"\n ");
    }
    fprintf(fp,"\n\n");
    fprintf(fp,"CELLS %d %d\n",n-1, 3*(n-1));
    for(int i = 0; i < n-1; i++)
    {
        fprintf(fp,"2 %d %d\n",i, i+1);
    }
    fprintf(fp,"\n\n");
    fprintf(fp,"CELL_TYPES %d\n", n-1);
    for(int i = 0; i < n-1; i++)
    {
        fprintf(fp,"3\n");
    }
    fclose(fp);
}

void fiq::vtkoutput(const char* str)
{
    double* pointer1 = this->f_inphase;
    double* pointer2 = this->f_quadrature;
    FILE* fp = fopen(str,"w");
    fprintf(fp,"# vtk DataFile Version 2.0\nVolume example\nASCII\nDATASET STRUCTURED_POINTS");
    fprintf(fp,"\nDIMENSIONS %d %d 1\nASPECT_RATIO %f %f 1\nORIGIN %.2f %.2f 0\n",
    Nx, Ny, dr, dz, Rmin, Zmin);
    fprintf(fp,"POINT_DATA %d\nSCALARS icompo float 1\nLOOKUP_TABLE default\n",
     Nx*Ny);

    for(int i = 0; i < Nx*Ny; i++)
    {
        int iy = i%Nx;
        int ix = i/Nx;

        fprintf(fp," %1.3f",pointer1[iy*Ny+ix]);
    }
    fprintf(fp,"\nSCALARS qcompo float 1\nLOOKUP_TABLE default\n");

    for(int i = 0; i < Nx*Ny; i++)
    {
        int iy = i%Nx;
        int ix = i/Nx;

        fprintf(fp," %1.3f",pointer2[iy*Ny+ix]);
    }
    fclose(fp);
}